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A  system  of  dynamical  equations  is  presented  that  allow  micro  air  vehicle  (MAV)  or- 
nithopter  designers  to  match  drive  motors  to  loads  produced  by  flexible  flapping  wing 
spars.  The  model  can  be  used  to  examine  the  coupled  system-level  behavior  of  brushed 
DC  motors,  gear  trains,  and  any  number  of  linkages  and  flexible  wing  spars.  A  Lagrangian 
approach  is  used  to  derive  the  governing  differential  equations  of  motion  for  a  class  of  or- 
nithopter  drive  systems.  Methods  used  to  determine  parametric  constants  contributing  to 
generalized  force  components,  which  cannot  be  derived  from  first  principles,  are  described. 

An  example  is  presented  where  simulation  results  are  compared  to  experimental  measure¬ 
ments.  The  results  show  that  the  differential  equations  correctly  predict  major  trends 
in  observed  motor  speed  and  wing  spar  structural  deformation  over  the  course  of  each 
wingbeat.  The  results  show  that  when  pairing  flight-weight  motors  and  wings,  significant 
variations  in  drive  motor  speed  occur  throughout  each  wingbeat.  It  is  shown  that  cou¬ 
pling  between  motor  speed,  wing  loads,  and  structural  flexibility  cause  the  aerodynamic 
forces  encountered  by  the  wing  spars  to  depart  from  those  predicted  by  rigid-spar  and 
constant-velocity-motor-based  kinematic  simulations. 

I.  Introduction 

A  large  class  of  micro-air  vehicle  (MAV)  ornithopter  flapping  mechanisms  consist  of  a  brushed  DC  motor 
and  gear  train  in  combination  with  linkage  elements,  flexible  wing  spars,  and  wing  surfaces.  The  aerodynamic 
and  inertial  loads  produced  by  the  wings,  and  the  inertial  loads  produced  by  gear  train  and  linkage  elements, 
can  significantly  affect  the  instantaneous  speed  of  the  motor  that  drives  them.  Improperly  matched  motors 
and  loads  can  result  in  repeated  starting  and  stopping  of  the  motor  armature,  which  limits  flapping  frequency 
and  reduces  system  performance.  Thus,  a  systematic  procedure  for  modeling  the  dynamics  of  this  class  of 
electro-mechanical  aeroelastic  systems  is  presented.  One  potential  use  for  such  a  model  is  to  serve  as  the 
basis  for  a  vehicle  design  tool  that  matches  drive  motors  to  loads  produced  by  flexible  wings.  A  major 
objective  of  this  exposition  is  to  present  a  low-order  dynamical  model,  of  adequate  fidelity,  to  predict  the 
characteristics  of  the  motion  of  flexible  wing  spars  that  are  driven  by  ornithopter  linkages  and  brushed  DC 
motors. 

There  are  three  principal  types  of  loads  that  a  drive  motor  must  overcome  in  an  ornithopter  application, 
viz.,  inertial,  frictional,  and  aerodynamic.  The  inertial  loads  arise  from  both  rigid  body  elements  and  flexible 
modes  associated  with  the  structural  dynamics  of  the  wings.  The  frictional  loads  arise  from  a  combination 
of  viscous  and  dry  friction,  where  the  former  is  due  to  light  lubricating  oils  on  the  gear  train  and  linkage 
elements.  Additionally,  the  back-electromotive  force,  or  back-EMF,  inherent  in  all  brushed  DC  motors 
produces  a  characteristic  that  from  the  point-of-view  of  the  motor  drive  torque,  behaves  like  viscous  friction. 
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The  Lagrangian  method  is  used  to  derive  a  system  of  governing  differential  equations  that  predict  the  coupled 
behavior  of  a  brushed  DC  motor  driving  a  gear  train  with  an  arbitrary  number  of  stages,  coupled  to  any 
number  of  linkages  driving  any  number  of  flexible  wing  spars.  Methods  for  obtaining  the  generalized  forces 
associated  with  each  degree-of-freedom  are  presented  for  the  non-conservative  forces  and  drive  elements 
within  the  system. 

Research  in  the  area  of  flapping  wing  flight  has  been  conducted  by  both  the  biological  and  engineering 
communities.  These  communities  have  leveraged  one  another’s  work  and  it  is  not  uncommon  for  biologists  to 
make  use  of  aerodynamic  theory  and  theoretical  mechanics  to  characterize  and  understand  the  mechanisms 
that  enable  flight  in  insects  or  birds. 1-4  Likewise,  the  engineering  community  has  taken  advantage  of  the 
work  of  biologists  to  both  devise  machines  that  can  approximate  the  motion  of  wings  observed  in  nature 
and  to  control  those  motions  in  such  a  way  as  to  enable  a  flapping  winged  aircraft  to  be  stabilized  and 
maneuvered  as  desired.  Shyy  et.  al.5  provides  an  overview  of  the  work  of  both  communities  and  provides  a 
foundation  for  understanding  the  physics  of  flapping  flight  at  low  Reynolds  numbers. 

The  current  state  of  technology  allows  one  to  design  and  construct  ornithopters  using  ad-hoc  methods 
to  match  loads  to  motors.  A  large  amount  of  the  work  reported  in  the  literature  describes  the  construction 
of  flying  prototypes.  For  instance,  the  flapping  wing  aircraft  known  as  DelFly6  was  developed  at  TU  Delft 
in  the  Netherlands  and  is  capable  of  controlled  forward  flight  and  hover.  A  similar  commercially  available 
aircraft,  known  as  the  WowWee  Dragonflya,  is  an  ornithopter  that  makes  use  of  four  flapping  wing  surfaces  to 
provide  propulsion  for  the  aircraft.  Researchers  at  the  Naval  Postgraduate  School  developed  an  aircraft  that 
replaces  conventional  propulsion  mechanisms'  with  flapping  wings.  Recently,  Gerdes  et.  al 8  exhaustively 
enumerated  the  flapping  wing  aerial  vehicles  that  have  had  at  least  one  successful  reported  flight.  In  a 
preliminary  study9  on  tailless  hover-capable  ornithopter  design,  it  was  found  that  such  aircraft  present 
major  challenges  associated  with  trim  and  control.  The  work  presented  here  shows  the  importance  of 
characterizing  interactions  between  loads  and  drive  motors,  and  the  effects  that  these  interactions  have 
upon  the  aerodynamic  forces.  It  is  important  for  control  law  designers  to  be  aware  of  the  presence  of  these 
interactions  and  their  impact  upon  the  cycle-averaged  behavior  of  the  aerodynamic  forces  and  moments. 

In  terms  of  designing  drive  mechanisms  for  flapping  wings,  it  is  desirable  to  minimize  the  number  of 
on-board  actuators  in  order  to  minimize  gross  take-off  weight.  Closed- loop  mechanisms  such  as  four-bar 
linkages  have  the  desirable  quality  of  creating  flapping  or  rocking  motion  using  a  single  drive  actuator. 
Extensive  work  has  been  reported  in  this  area.  Researchers  have  used  synthesis-based  or  optimization-based 
methods  to  design  mechanisms  that  match  the  wing  motion  produced  by  many  flying  biological  species  in 
a  kinematic  sense.  For  instance,  the  series  of  designs  by  Agrawal,  et.al. 10,11  attempted  to  optimize  various 
mechanism  designs,  based  upon  kinematic  or  quasi-static  models  of  the  mechanism  motion,  in  order  to  match 
the  motion  generated  by  some  of  the  representative  flying  species.  Similarly,  Rahmat  et.  al 12  synthesized  and 
fabricated  a  spherical  four-bar  mechanism  to  mimic  the  figure-8  wing  stroke  motion  observed  in  a  number 
of  natural  flyers.  While  they  showed  simulation  and  experimental  results,13  their  proof-of-concept  prototype 
only  qualitatively  showed  the  feasibility  of  such  design.  Raney  and  Slominski14  developed  a  ground-test 
apparatus  that  could  vary  the  wing  kinematics  produced  by  a  machine  that  approximated  the  wing  motion 
observed  in  hummingbirds.  The  experiment  focused  on  the  control  of  the  wingtip  trajectories  by  varying 
the  amplitude  and  phase  of  waveforms  applied  to  two  electrodynamic  actuators.  Chung  et.al.15  investigated 
the  use  of  central  pattern  generators  to  control  the  three  dimensional  motion  of  the  wings  of  a  10  degree-of- 
freedom  robotic  bat  in  a  wind  tunnel  experiment.  Other  researchers16-18  have  investigated  the  problem  of 
optimizing  wingbeat  kinematics  for  power  or  aerodynamic  efficiency. 

In  the  above  efforts,  however,  the  interactions  between  the  drive  elements  and  wing  flexibility  were  not 
mathematically  modeled  when  optimizing  mechanisms  for  matching  the  observed  motion  of  natural  flyers. 
In  the  present  exposition,  it  is  shown  that  the  idealized  kinematic  behavior  of  a  linkage  and  rigid  wing  spar 
driven  by  a  constant  velocity  DC  motor  is  quite  different  from  that  produced  by  a  system  where  a  torque- 
limited  flight- weight  motor  drives  flexible  wing  spars.  As  will  be  shown  in  the  present  work,  the  ability  to 
prescribe  wingbeat  kinematics  in  flight-weight  systems  will  not  normally  be  possible  because  of  torque  limits 
and  dynamic  interactions  between  drive  and  wing  components.  One  of  the  conclusions  of  the  present  study 
is  that  the  entire  system,  from  drive  motor  to  wing,  must  be  integrated  into  the  model  when  attempting  to 
match  the  observed  motion  of  biological  systems  with  a  flight-weight  mechanical  system. 

In  terms  of  overall  flight  control,  Deng  et.al.19,20  designed,  modeled,  simulated,  and  produced  a  cycle- 
averaged  control  design  for  a  conceptual  biomimetic  flapping  wing  robotic  aircraft.  Khan  and  Agrawal21 

ahttp:  /  /  www.wowwee.com /en/ products  /  toys /flight /flytech  /  dragonfly 
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used  a  time-averaged  control  approach  in  conjunction  with  a  differential  flatness  based  control  scheme  to 
design  a  longitudinal  flight  controller  for  a  flapping-wing  MAV.  DeLeo  and  Deng22  presented  a  four-wing 
aircraft  design  based  on  a  dragonfly  that  featured  a  single  rotary  motor  that  drove  the  wings  in  the  stroke- 
plane,  a  tunable  phase  difference  between  the  forewings  and  hindwings  using  a  variable  length  coupler  link, 
and  passive  wing  rotation  about  the  spar.  In  the  above  control  related  studies,  the  dynamic  interactions 
between  the  drive  components  and  the  flexible  and  aerodynamic  loads  were  not  considered.  Because  such 
interactions  affect  the  aerodynamic  forces  and  moments  produced  by  the  wings,  they  will  also  affect  one’s 
ability  to  control  the  aircraft. 

Finally  it  is  worth  noting  that  there  is  line  of  research  pursued  by  Wood23  et.al.  that  has  focused  on 
the  development  of  piezoelectrically  powered  micro  mechanical  insects.  The  first  takeoff  of  an  insect-scale 
flapping  wing  MAV  was  achieved  by  an  aircraft  called  RoboFly  that  was  developed  at  Harvard  University 
by  Wood  et.al.24  Avadhanula,  et.al.25  modeled  a  flapping  wing  piezo-electrically  actuated  electromechanical 
subsystem  system  using  a  Lagrangian  approach.  While,  piezoelectrically  actuated  flapping  wing  aircraft  are 
an  important  class  of  vehicles  with  a  promising  future,  the  work  presented  here  is  focused  on  the  funda¬ 
mentally  different  dynamic  interactions  that  occur  between  DC  motor  powered  gear-trains,  linkage-based 
mechanisms,  and  flexible  wing  spars. 

The  remainder  of  the  paper  is  structured  as  follows:  Section  II  presents  a  description  of  the  components 
that  comprise  the  class  of  systems  under  consideration.  The  Lagrangian  formalism  enables  the  analysis  of  the 
interactions  between  each  subsystem,  in  terms  of  energy,  and  is  presented  in  Section  III.  Section  IV  evaluates 
the  equation-of-motion  of  the  system  using  Lagrange’s  equation.  Taking  advantage  of  the  additive  properties 
of  derivatives,  the  equations-of-motion  are  derived  in  a  modular  fashion.  In  Section  V,  nonconservative 
generalized  forces,  such  as  drive-motor  torque,  modal  damping,  frictional,  and  aerodynamic  forces  are  derived 
using  the  principle  of  virtual  work.  The  information  required  to  evaluate  these  terms  are  supported,  in  some 
cases,  by  experimental  measurements.  Finally,  in  Section  VI,  simulation  results  are  critically  compared  to 
experimental  results  generated  by  the  system  described  in  Section  II. 

II.  Description  of  a  Bench-Test  Article 

The  framework  developed  in  this  manuscript  is  directly  applicable  to  a  large  class  of  MAV  ornithopter 
designs;  however,  the  simulation  and  experimental  results  that  are  presented,  apply  to  a  specific  test  article 
that  consists  of  components  of  a  prototype  aircraft  considered  by  Doman  et.al.9  In  order  to  introduce  the 
reader  to  the  class  of  systems  under  consideration,  a  specific  bench  test  article  will  now  be  described.  The 
bench-test  article  and  some  of  the  supporting  measurement  apparatus  is  shown  in  Fig.  1.  For  the  purpose  of 
the  analysis  in  this  paper,  the  system  is  considered  to  be  rigidly  mounted  to  the  bench.  The  system  consists 
of  a  coreless  brushed  DC  motor  driving  a  set  of  compound  gears,  which  in-turn  drives  a  four-bar  mechanism 
that  imparts  rocking  motion  to  a  flexible  wing  spar. 


(a)  Bench-test  article  (b)  Enlarged  view  of  driving  mechanism  (diameter  of  large 

gear  is  2.49  cm) 


Figure  1.  Experimental  Setup 
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The  fuselage  frame  was  produced  from  a  CAD  model  using  a  3-D  printer.  The  gears  came  from  a 
commercial  source.  The  linkage  elements  were  custom  machined  from  Delrin.  The  linkage  pins  were  made 
from  commercially  available  1/16  inch  (.158  cm)  aluminum  rivets.  The  spar  was  made  from  commercially 
available  0.037  inch  (.09398)  carbon  fiber  rod  cut  to  13cm  in  length.  The  spar  was  attached  to  the  rocker 
arm  of  the  linkage  using  a  Kevlar  wrap  that  was  saturated  with  cyanoacrylate  adhesive.  The  fuselage  frame 
also  incorporated  a  custom-made  16  slot  optical  encoder.  The  encoder  meshed  with  the  linkage  crank  gear 
to  measure  the  overall  system’s  angular  velocity.  Because  of  the  choice  of  gear  ratios,  the  encoder  rotated 
at  the  same  velocity  as  the  stage  1  gear.  The  motor  power  was  supplied  by  a  regulated  voltage  source.  A 
high  speed  camera  system  was  used  to  capture  the  motion  of  the  system  at  2000  frames  per  second  and 
the  images  were  used  to  critically  compare  the  motion  of  the  simulation  model  with  that  observed  in  the 
experiments. 

The  development  that  follows  applies  to  a  general  class  of  ornithopter  flapping  mechanisms.  The  bench- 
test  mechanism  described  above  was  used  to  experimentally  validate  the  equations  of  motion  for  a  special 
case. 


III.  Modeling  and  Lagrangian  Formulation  of  Each  Subsystem 

The  overall  system  consists  of  the  following  components:  a  motor,  gear-train,  four-bar  linkage,  and 
a  flexible  wing  spar.  For  the  purpose  of  this  paper,  the  fuselage  is  assumed  to  be  rigidly  mounted  to  a 
bench,  thus  a  fuselage-fixed  coordinate  frame  is  taken  to  be  an  inertial  frame.  The  motor  and  the  linkage 
line-of-centers  are  also  taken  to  be  fixed  with  respect  to  the  fuselage  frame. 

In  what  follows,  the  kinetic  energy  (T)  and/or  potential  energy  (V)  of  each  subsystem  is  considered.  The 
Lagrangian  of  a  mechanical  system  is  defined  as: 

C  =  T-V  (1) 

In  a  subsequent  section,  the  energy  of  all  of  the  subsystems  are  summed  together,  and  the  equations-of-motion 
of  the  overall  system  are  determined  using  Lagrange’s  equation. 

A.  DC  Motor 

First,  a  model  of  a  brushed  DC  motor  is  considered.  According  to  Ogata,26  the  torque  produced  by  a  motor 
is  proportional  to  the  current  flowing  through  the  windings: 


ts  =  KTia  (2) 

where  Kt  is  the  motor  torque  constant,  and  ia  is  the  current.  The  back-EMF  produced  by  the  motor  as  a 
result  of  its  angular  velocity  can  be  described  as: 

e6  =  Kb9  (3) 

where  Kb  is  the  back-EMF  constant,  9  is  the  angular  position  of  the  motor  armature,  and  9  is  its  angular 
velocity.  The  differential  equation  governing  the  dynamics  of  the  armature  circuit  is  given  by: 

dia 

La  H-  Llaia  H-  65  —  ca  (4) 

at 

where  La  is  the  inductance  of  the  motor,  Ra  is  the  resistance  of  the  armature,  and  ea  is  the  voltage  applied 
to  the  armature.  According  to  Ogata,26  motor  inductance  is  typically  very  small;  thus,  La  ss  0.  Using  the 
small  inductance  approximation  and  solving  for  ia  from  Eq.  (4)  yields: 


la  = 


&CL 


Rn 


(5) 


Substituting  Eqs.  (5)  and  (3)  into  Eq.  (2)  allows  one  to  solve  for  the  net  torque  produced  by  the  motor, 


Kt  KTKb  ■ 

tb  =  —ea - —9 

■Tt  a  -tv  n, 


(6) 
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Note  that  the  torque  is  dependent  on  the  input  voltage  e0.  The  back-EMF  and  drive  torques  form  a  portion 
of  the  total  generalized  force  associated  with  the  rigid  body  degree-of-freedom  6. 

The  rotation  of  the  motor  armature  contributes  kinetic  energy  to  the  mechanical  system  according  to: 

Tm  =  \jae2  (7) 


where  Ja  is  the  mass  moment  of  inertia  of  the  armature  and  spindle  about  the  axis-of-rotation.  Since  the 
motor  does  not  contribute  to  the  potential  energy  of  the  system,  the  Lagrangian  of  the  DC  motor  subsystem 
is  given  by: 


r  —  T  —  -  T 

—  J-m  —  2  ,Ja 


(8) 


B.  Gear  Train 

Many  ornithopter  designs  use  compound  gear  trains  to  reduce  the  speed  of  the  motor  to  achieve  a  target 
wing  flapping  frequency.  The  overall  speed  reduction  and  torque  amplification  can  be  tuned,  within  limits, 
to  match  the  torque  of  the  motor  to  the  loads  produced  by  the  wing  and  the  mechanism. 


Spur  Gear  n 
or  Crank  Gear 


\e,  \  =  r  \e 


Kinetic  Energy  of 
Each  Gear 


Tg 

ong  'j  on 


i  'b-1 

T  =-T[J  Gfd2 

Sn- 1  9  X  X  Sk-l  k 

z  k= 1 


r  =  —j  G2G292 

S  2  2  82  2  i 


T ’  =-J  G292 

gi  2  Sl 


Motor  T  =  —J  61 

go  2  g° 


Figure  2.  Compound  Gear  Train  Transmission 


Figure  2  shows  a  compound  gear  train  that  consists  of  ng  stages.  The  stage  0  gear  in  the  train  is  a  pinion 
fixed  to  the  motor  armature.  The  armature  pinion,  which  has  po  teeth,  meshes  with  periphery  of  a  first 
stage  compound  gear  that  has  t\  teeth.  At  each  stage,  the  pinion  of  each  compound  gear  meshes  with  the 
spur  gear  of  the  subsequent  stage  to  incrementally  increase  mechanical  advantage  and  reduce  speed  until  the 
final  crank  gear  stage  ng  is  reached.  The  crank  gear  stage  ng  then  forms  the  crank  for  the  four-bar  linkage 
portion  of  the  drive-train.  The  total  speed  reduction  from  the  pinion  to  the  crank  is  given  by: 


e, 


Ifa 


vfc= o 


(9) 
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where  9i  is  the  crank  angular  velocity  of  the  four-bar  linkage,  ng  is  the  number  of  gears  in  the  gear  train, 
Gk  is  the  gear  ratio  of  the  kth  stage,  which  for  a  given  gear  pitch,  is  given  by: 


Gk  = 


1, 

Pk- 1 

tk  ’ 


k  =  0 

k=  1,2,--- 


(10) 


where  pu-i  is  the  number  of  teeth  on  the  pinion  of  compound  gear  k  —  1,  and  tk  is  the  number  of  teeth  on 
the  spur  gear  k.  For  convenience,  the  overall  gear  ratio  from  the  motor  to  the  crank  gear  is  defined  as: 


ng 

r  =  n  Gk 


k=0 


(ii) 


Excluding  the  crank  gear,  the  total  kinetic  energy  and  the  Lagrangian  of  the  ng  —  1  gears  and  the  pinion 
attached  to  the  motor  spindle,  is  given  by: 


ng  —  1  fc 


E 

k— 0  j—0 


(12) 


where  Jgo  is  the  mass  moment  of  inertia  of  the  pinion,  and  Jgk  is  the  mass  moment  of  inertia  of  each 
compound  gear  in  the  gear  train.  Note  that  the  crank  gear  energy  is  not  included  in  Tg  because  it  is 
included  in  the  energy  of  the  linkage  that  is  derived  in  the  subsequent  section.  Since  there  is  no  potential 
energy  contribution  from  the  gear-train,  the  Lagrangian  of  the  gear  train  subsystem  is  simply  Cg  =  Tg: 

For  the  system  introduced  in  Section  II,  ng  =  2.  Hence: 

£g  =  Tg  =  ^(JgO  +  JglGl)@2  (13) 

It  is  worthwhile  to  note  that  the  kinetic  energy  of  other  rotating  components,  such  as  encoder  wheels  used 
to  measure  position  or  speed,  can  be  included  by  lumping  the  mass  moment  of  inertia  of  these  components 
with  the  mass  moment  of  inertia  of  any  gears  that  rotate  at  the  same  speed.  For  example,  the  encoder  wheel 
used  in  the  experiments  performed  under  this  work,  moves  at  the  same  speed  as  the  first  stage  compound 
gear.  Thus,  the  effective  mass  moment  of  inertia  can  be  computed  from  the  following  sum: 


Jgl  —  Jgl  +  Je  +  J ep  (14) 

where  Je  and  Jep  represent  the  mass  moments  of  inertia  of  the  encoder  wheel  and  the  pinion  attached  to  the 
encoder  wheel  respectively. 


C.  Four-Bar  Linkage 

Figure  3  defines  the  elements  of  a  four-bar  linkage  for  the  class  of  systems  under  consideration.  The  physical 
structure  of  the  crank  gear,  viz.,  gear  stage  ng,  forms  the  crank  of  the  four-bar  linkage.  A  pin  of  mass  mp 
connects  the  crank  gear  to  a  coupler  bar  of  mass  m2  and  mass  moment  of  inertia  about  its  own  center-of- 
mass  of  I2  -  The  coupler  bar  is  connected  to  a  rocker  arm  of  mass  m3  and  mass  moment  of  inertia  J3  by  a 
pin  of  mass  rnp.  The  line  that  connects  the  center-of-rotation  of  the  crank  gear  and  the  pivot  point  of  the 
rocker-arm  is  called  the  line-of-centers,  whose  length  is  denoted  as  Iq. 

Due  to  the  closed-loop  kinematic  constraints  formed  within  the  linkage,  the  time  evolution  of  the  gen¬ 
eralized  coordinates  of  the  crank  gear  angle  9i,  coupler  angle  cq,  and  rocker  angle  (j)i  are  highly  coupled. 
However,  it  is  well-known  that  a  four-bar  linkage  is  a  single-degree-of-freedom  system,  permitting  one  to 
express  cq  and  </>;  in  terms  of  9\.  Hence,  the  objective  of  this  section  is  to  establish  the  total  kinetic  energy 
of  all  the  linkage  elements  in  terms  of  the  motor  angular  velocity  9. 


1.  Four-Bar  Linkage  Kinematics 

To  solve  the  kinematic  problem  of  the  four-bar  linkage  under  consideration,  the  loop-closure  constraint 
equation  is  written  as: 

li  cos  9i  +  Z2coscq  =  la  +  l3  cos  cj)i  (15) 

Zisinf?/  +  Z2sincq  =  /3sin0;  (16) 
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Figure  3.  Four  Bar  Linkage 


Solving  for  the  terms  containing  a\  yields: 

1-2  cos  ai  =  l0  +  h  cos  <j>i  —  l\  cos  9i 
I2  sin  ai  =  l3  sin  (f>i  —  l\  sin  9i 

Eliminating  oti  by  squaring  and  summing  Eqs.  (17)  and  (18)  yields  the  Freudenstein  Equation: 

I<i(9i)  sin  fa  +  K2(9i)  cos  4n  +  K3{9i)  =  0 


where, 


Ki{9i)  =  —2l\l3  sin  9i 
K2(9i)  =  2 l3(l0  -  h  cos 9i) 

K3(9i)  =  l2a  +  li  —  l\  +  Zg  —  2loh  cos  9i 


Letting  r  =  tan(</>;/2),  and  making  use  of  double-angle  trigonometric  identities  yields: 


sin  (j>i 
cos  <fn 


2  r 

1  +  T2 
1-T2 

1  +  r2 


(17) 

(18) 


(19) 


(20) 

(21) 

(22) 


(23) 

(24) 


Substituting  Eqs.  (23)  and  (24)  into  Eq.  (19)  yields  a  quadratic  equation  in  r  for  a  given  crank  angle  9f. 


(it's  -  K2)t2  +  2K\t  +  (K3  +  I<2)  =  0 


(25) 


from  which  one  can  readily  solve  for  the  rocker-arm  angle  <f>i  via: 


<t>i  =  2  •  arctan2 


'-K!±  y/Kf-iq+K*' 


K3  -  Ko 


(26) 


where  arctan2  denotes  the  4-quadrant  arc-tangent  function  that  is  dependent  upon  the  signs  of  the  numerator 
and  denominator.  The  coupler  angle  oti  can  then  be  solved  by  dividing  Eq.  (17)  by  Eq.  (18): 


oti  =  arctan2 


(27) 


Hence,  ai(9i)  and  <f>i(8i)  can  effectively  be  written  as  functions  of  9/. 
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The  relationships  between  the  crank  angular  velocity  9i,  and  the  dependent  angular  velocities,  a i  and 
<j)i,  can  also  be  determined.  Differentiating  the  loop  closure  equations  in  Eqs.  (15)  and  (16)  with  respect  to 
time  yields: 

— l\6i  sin  61  —  I2011  sin  cq  +  l3(f>i  sin  <j>i  =  0  (28) 

l\0i  cos  9i  +  l2dn  cos  cq  —  l3<fii  cos  <j>i  =  0  (29) 


Equations  (28)  and  (29)  can  be  conveniently  written  as  a  linear  algebra  problem  that  can  be  solved  for  cq 
and  <t>i  via  Cramer’s  rule  yielding: 


h  sin  (0;  -  9i)  • 

ai  =  - Jr0* 

h  sin  (cq  -  4n ) 

■  _  h  sin  (cq  -  9i)  ■ 

1  l3  sin  (cq  -  00  1 

Since  9i  =  T9,  Eqs.  (30)  and  (31)  can  be  written  as: 

ai  =  S!(9)6=^e 
j>i  =  32(9)9  =  ^ 9 

where, 

c  (Q\  Arli  sin  (0;  -  T9) 
$1  (9)  —  1  T^l - TT 

1 2  sm  ( cq  -  0/ ) 
c  ta\  a  sin  (at  -  T9) 

J2(V)  =  1  7 - TT 

<3  sm(tt,  -  0/) 


(30) 

(31) 

(32) 

(33) 

(34) 

(35) 


Therefore,  the  linkage  element  velocities  A;,  0;,  and  9i  can  be  written  in  terms  of  a  single  motion  variable 
and  its  time  derivative,  viz.  the  the  motor  angular  position  9  and  angular  velocity  9. 


2.  Lagrangian  of  Four-Bar  Linkage 

The  kinetic  energy  of  the  linkage  elements  is  given  by: 


T'  =  \ 


willed  ||2  +  h9f  +  to2||uC2||2  +  I2af  +  m3||uC3||2  +  J302  +  mp(||uPl  ||2  +  ||uP2||2) 


(36) 


where  vCl ,  vC2  and  vC3  are  the  linear  velocities  of  the  centers-of-mass  of  the  crank,  coupler,  and  rocker 
respectively;  rnp  is  the  mass  of  the  linkage  pins,  and  vPl  and  vP2  are  the  velocities  of  the  two  pins.  The 
center-of-mass  of  the  crank  is  constant  because  it  is  formed  from  the  final  gear,  which  is  axis-symmetric, 
thus  vCl  =  0  and  I3  =  Jgrlg .  Expanding  Eq.  (36)  yields: 


l  r 


T'=2 


Ii&f  +  m,2(l\9i  +  l2C2af  +  2l\lC2  cos (9i  —  cq)cq#;)  +  I2A2  +  (^3!^  +  /3)02  +  +  Z302)  |  (37) 


Collecting  terms  on  the  products  of  the  angular  velocities  and  substituting  Eqs.  (9)  and  (11)  into  Eq.  (37) 


yields: 

Ti  —  -[Jij02  +  J2loq  +  J3i02  +  2P\C\ai9] 

(38) 

where, 

Jq  =  T2  [/i  +  m2l\  + 

(39) 

■J‘2i  =  JTl2^c2  + 

(40) 

J$i  =  ni3 l23  +  mpl3  +  I3 

(41) 

P\  =  m2hlc2  T 

(42) 

Cl  =  cos (9i  -  cq) 

(43) 
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Note  that  the  linkage  operates  in  a  plane  that  is  perpendicular  to  the  gravitational  vector,  thus,  there  is 
no  change  in  potential  energy  of  the  linkage  elements  as  the  crank  rotates.  Hence,  only  kinetic  energy 
contributes  to  the  total  Lagrangian  associated  with  the  four-bar  linkage  subsystem.  Substituting  Eqs.  (30) 
and  (31)  into  Eq.  (38)  yields  the  complete  expression  for  the  kinetic  energy  and  Lagrangian  of  the  linkage 
in  terms  of  the  motor  angular  velocity  9: 

Ci  =  Tt=±  [Ju  +  J2l  Sf  +  J3l  Si  +  2P1C1S1}92  (44) 


D.  Flexible  Spar 

Typical  MAV  ornithopters  use  carbon  fiber  components  for  the  wing  structure  and  thin  films,  such  as  Mylar, 
for  the  wing  skin.  In  general,  the  interactions  of  the  motor,  gears,  linkage,  and  wing  are  of  interest.  In  the 
analysis  that  follows,  the  structural  dynamics  of  a  wing  spar  are  considered,  since  the  leading  edge  wing  spar 
provides  the  primary  resistance  to  wing  bending  in  the  stroke  plane.  Figure  4  shows  a  diagram  of  a  four  bar 
linkage  with  a  flexible  wing  spar  clamped  to  the  end  of  a  rocker-arm. 


Figure  4.  Flexible  spar  driven  by  rocker  arm. 


An  inertial  frame,  represented  by  the  unit  vectors  hi,fi2,  is  fixed  at  the  center-of-rotation  of  the  rocker- 
arm.  A  coordinate  frame  that  rotates  with  the  rocker  arm  is  denoted  as  si,  S2-  A  vector  from  the  origin  of 
frame  si  —  S2  to  the  point  at  which  the  spar  is  clamped  to  the  rocker  arm  is  denoted  as  r>.  A  vector  from 
the  clamped  root  of  the  spar  to  an  arbitrary  point  on  the  undeformed  spar  is  denoted  as  ru.  A  vector  from 
a  point  on  the  undeformed  spar  to  a  point  on  the  deformed  spar  is  denoted  as  w(ru,t )  where  ru  =  |ru|  . 
The  deformation  w(ru,t )  is  taken  to  be  perpendicular  to  the  undeformed  spar  and  therefore  occurs  only  in 
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the  §2  direction.  When  viewed  from  the  rocker-arm  frame,  the  flexible  spar  appears  as  a  cantilever  beamb. 
The  motion  of  the  rocker-arm  will  cause  the  spar  to  deform  due  to  bending.  Euler-Bernoulli  beam  theory  is 
used  in  the  following  analysis  and  the  spar  is  taken  to  have  uniform  cross  section  and  Young’s  modulus.  The 
assumed  modes  method27  is  used  to  describe  the  deformation  of  the  spar  in  the  stroke  plane  relative  to  the 
si  axis.  The  deformation  is  taken  to  be  described  by  an  infinite  sum  of  the  products  of  spatially  dependent 
mode  shape  functions  &i(ru)  and  temporal  modal  coordinates  i.e.: 


w(ru,t)  =  ’y]$i(ru)rii(t) 

i= 1 

where,  for  clamped-free  boundary  conditions,  the  ith  mode  shape  function  is  given  by: 

$i(r„)  =  Ai  [aj(sin/3,ru  -  sinh  /^ry)  +  &j(cos/?;ru  -  cosh^ru)] 


where, 


a,;  =  (sin  fiiL  —  sinh  piL) 
bi  =  (cos  PiL  +  cosh  PiL) 


(45) 


(46) 


(47) 

(48) 


and  Pi  is  the  ith  solution  to  the  frequency  equation: 

cos  piL  cosh  piL  —  —1  (49) 


and  Ai  is  an 
its  length: 


arbitrary  constant  with  units  of  length0.  For  a  cantilever  beam  with  uniform  properties  along 


Pl  = 


ufrh 

E„h 


(50) 


where  uy  is  the  natural  frequency  of  the  ith  bending  mode,  m  is  the  mass  per  unit  length  of  the  beam,  Ea  is 
the  Young’s  modulus  of  the  spar,  and  Is  is  the  area  moment  of  inertia  of  the  cross  section  of  the  spar  about 
the  neutral  axis. 

In  the  analysis  that  follows,  the  first  and  second  bending  modes  are  retained  in  the  response11.  Hence, 


w(ru,t )  =  $i(ru)rji(t)  +  $2(ru)V2(t) 
Solving  the  frequency  equation  for  the  first  two  bending  modes  yields: 


uq  =  1.875" 


u>2  =  4.694" 


ESIS 

rhLA 

ESIS 

rhLA 


where  the  constants  are  accurate  to  4  significant  figures.  It  is  also  useful  to  note  that: 


01 

02 


1.875 

L 

4.694 

L 


(51) 

(52) 

(53) 

(54) 

(55) 


1.  Flexible  Spar  Kinetic  Energy 

For  convenience,  define  r  =  ry  +ru.  Note  that  rr  extends  along  the  rigid  rocker  arm  from  the  pivot  point  to 
the  clamped  end  of  the  spar.  From  Fig.  4,  one  may  write  the  position  vector  from  the  origin  of  the  rocker 
arm  to  a  point  on  the  deformed  beam  as: 

rp  =  [r  cos  pi  —  w  sin  pi]  h\  +  [r  sin  pi  +  w  cos  pi]  fi2  (56) 

“Note  that  the  rocker  arm  angle  <pi  remains  a  function  of  motor  position  8 

cIn  the  literature,  the  arbitrary  constant  associated  with  each  mode  is  often  used  to  mass  normalize  the  mode  shapes 
dFor  the  case  considered  in  the  results  section,  it  was  found  that  retaining  only  the  first  bending  mode  yielded  excellent 
results;  nevertheless,  a  more  general  two-mode  case  has  been  considered  in  the  analytical  development 
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Substituting  Eq.  (51)  yields: 


rp  =  [r  cos  <f>i  -  ($ir?i  +  $2?72)  sin  fa]  hi  +  [r  sin  fa  +  ($iT?i  +  $2*72)  cos  fa]n2 


(57) 


Differentiating  Eq.  (57)  with  respect  to  time  yields: 

rp  =  \-rfa  sin  fa  -  (<J>i??i  +  §2^2)  fa  cos  fa  -  ($i?ji  +  $2772)  sin  c 


ni  + 


rfa  cos  (j>i  -  ($i??i  +  $2?7 2)  fa  sin  ^  +  ($1771  +  $2?)2)  cos  < 


n2 


The  differential  kinetic  energy  of  an  elemental  beam  segment  on  the  deformable  spar  is  given  by: 


dTs  =  -rh[fp  ■  fp]dr 


or,  with  some  mathematical  manipulation 
1 


dTs  =  -?h 


r  +  ($1771  +  $2772)  </>f  +  ($i?7i  +  $2?72)  +  2r^($i?7i  +  $2772) 


dr 


Integrating  Eq.  (60)  from  the  root  to  the  tip  of  the  spar  yields: 


rrr+L  r 


Ts  =  2m 


f2  4>i  +  ($l?7l  +  $2772  )24>l  +  (^1771  +  ^27?2)2  +  %rfa  ($  1?71  +  <h2?72) 


\  2  12 


dr 


(58) 


(59) 


(60) 


(61) 


or 


rL  r 


Ts  =  -rh 
2 


(rr  +  ru)2j>f  +  (<hi?7i  +  <E>2772)2<^>f  +  (^it)!  +  <f>2?72)2  +  2(rr  +  ru)fa{^ifji  +  <&2?72)  dru  (62) 


Recall  that  the  shape  functions,  4>i,  are  functions  of  ru  only,  and  that  <pi,  771  and  772  are  functions  of  time 
and  do  not  vary  with  ru.  Carrying  out  the  integration  of  Eq.  (62)  yields: 


rr  1  ~ 

Ts  =  2m 


where: 


Cifa  +  c2{fa  ?7i  +  T7i )  +  2c40;?7i  +  2c6fair]2fa  +  771772)  +  c7(</>z  7?2  +  t)2)  +  2c8<(>;772 


7^  L6 

ci  =  /  (?y  +  ru)2dru  =  —  +  L27>  +  Lr2 

Jo  ^ 

c2  =  /  d>2(iru  =  9.226L 
Jo 

c4  =  /  (r>  +  rtl)<l>i(i7’11  =  — L(1.728L  +  2.378rr) 
Jo 


Ce  —  dru  =  0 

Jo 

c7=  [  fa2dru  =  2984.480L 

70 
rL 


Cs  =  f  (? >  +  ru)$2dru  =  — L(4.959L  +  23.707rr) 
Jo 


(63) 

(64) 

(65) 

(66) 

(67) 

(68) 

(69) 


Note  that  the  decimal  constants  in  c2  and  c7  have  units  of  length  squared,  while  the  decimal  constants  in 
c4  and  cs  have  units  of  length,  and  cq  is  0  since  the  mode  shapes  are  orthogonal.  Substituting  Eq.  (35)  and 
Eq.  (67)  into  Eq.  (63)  yields: 


rr  1  . 

Ts  =  -m 


CiS'le2  +  c2{S292ril  +  ?)2)  +  2cAS29f]i  +  c7(S292r j2  +  rfe)  +  2c8S29fi2 


(70) 


The  kinetic  energy  of  the  spar  is  now  written  in  terms  of  the  three  motion  variables  of  interest,  viz.,  9,  771, 
772,  and  their  derivatives. 
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2.  Flexible  Spar  Strain  Energy 

By  definition,  the  strain  energy  due  to  bending  of  an  Euler-Bernoulli  beam  is  given  by28  : 

VS=^ESISJ  w"(ru,  t)2dru  (71) 

Retaining  only  the  first  two  bending  modes  and  setting  the  arbitrary  constant,  A,;,  of  Eq.  (46)  to  unity,  the 
second  spatial  derivative  of  the  deformation  becomes: 


w"(ru,t)  =  Pi  [— ai(cos/?iru  +  cosh/?ir„)  -  &i(sin/3ir„  +  sinh/3ir„)]  ?yi+ 
Pl  [-a2(cos/32r„/  +  cosh p2ru)  -  b2(sinp2ru  +  sinh (32ru)]  q2 


(72) 


Squaring  Eq.  (72)  and  evaluating  Eq.  (71)  yields  the  strain  energy  of  the  spar  due  to  bending  deformation 
from  the  first  two  bending  modes: 

vs  =  7^EsIs  [c5r)l  +  Ci0??2 ]  (73) 

where, 


c  5 


ClO 


114.082 

L3 


1.449  x  106 

Z3 


(74) 

(75) 


Integrals  of  the  products  of  the  two  mode  shapes  are  zero  because  of  orthogonality.  The  decimal  constants 
in  C5  and  C10  have  units  of  length  squared.  The  contribution  of  the  flexible  spar  subsystem  to  the  Lagrangian 
is: 


ES  =  TS-  14 


(76) 


or 


Cs 


c\S\d2  +  c2(Sld2r]l  +  i)l)  +  2caS26pi  +  ct^S'IO2^  +  r]2)  +  2csS29f]2 


ESIS  [c5?y2  +  ci0?72]  (77) 


IV.  Equations-of-Motion 


In  Section  III,  the  Lagrangian  contribution  of  each  subsystem  was  evaluated.  The  full  equations-of-motion 
can  now  be  derived  using  Lagrange’s  equation: 


d_  fdC\ 

dt  \dqi ) 


(78) 


where  qi  represents  the  fh  generalized  displacement  coordinate,  and  Qi  the  corresponding  generalized  forces 
derived  from  the  principle  of  virtual  work. 

For  the  system  under  consideration,  there  are  3  generalized  coordinates,  namely,  0 ,  r/i,  and  ?y2-  Hence, 
three  coupled  ordinary  differential  equation-of-motions  are  derived.  The  Lagrangian  of  the  entire  system  can 
be  obtained  by  summing  Eqs.  (8),  (12),  (44),  and  (77): 


E  —  Em  +  Eg  +  Ei  +  Es  (79) 

While  this  paper  studies  a  system  where  a  DC  motor  drives  a  single  gear  train  and  four-bar  linkage 
assembled  with  a  flexible  spar  attached  to  the  rocker  arm,  it  is  straight-forward  to  extend  the  formulation 
to  account  for  systems  with  additional  elements.  For  instance,  it  is  common  for  ornithopters  to  use  a 
single  motor  to  drive  two  linkage-spar  assemblies,  as  presented  by  Doman  and  Regisford,9  or  a  single  linkage 
assembly  to  drive  two  spars.  For  such  cases,  let  rq  be  the  number  of  linkage  assemblies  and  ns  be  the  number 
of  spars  attached  to  the  linkages.  Assuming  that  only  identical  copies  of  such  elements  are  added  to  the 
system,  the  more  general  Lagrangian  can  be  written  as: 

E  =  Em  +  Eg  +  niEi  +  nsEs  (80) 
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The  remainder  of  this  section  evaluates  the  inertial  and  Coriolis/centrifugal  terms  associated  with  the 
left-hand  side  of  Lagrange’s  equation  in  Eq.  (78).  In  the  interest  of  establishing  a  general  set  of  equations 
in  which  components  may  be  added  or  subtracted  from  the  system,  we  note  that  the  derivative  of  a  sum  is 
equal  to  the  sum  of  the  derivatives,  thus  Lagrange’s  equation  can  be  evaluated  for  each  component,  and  the 
results  can  be  added  together  to  create  the  full  equations  of  motion.  One  can  derive  the  more  general  cases 
using  the  Lagrangian  shown  in  Eq.  (80). 


A.  DC  Motor 

Referring  to  Eq.  (8),  it  can  be  seen  that  the  Lagrangian  associated  with  the  motor  is  a  function  of  the  single 
motion  variable  9.  Hence,  it  is  straight-forward  to  show  that  the  contribution  of  the  motor  to  the  equation 
of  motion  associated  with  the  rigid  body  degree-of-freedom  is  given  by: 


d_  (dCm\ 
dt  V  99  ) 


=  JJ 


The  motor  also  contributes  to  the  nonconservative  input  torque: 


Q 


me 


re 


(81) 


(82) 


Recall  that  tq  includes  the  effects  of  the  voltage  applied  to  the  armature  and  the  back  EMF  damping  effect, 
which  is  a  function  of  the  motion  variable  9. 


B.  Gear  Train 


The  Lagrangian  expression  for  the  gear  train  in  Eq.  (12),  is  a  function  of  9  only.  Hence,  the  contribution  of 
the  gear  train  to  the  equation  of  motion  of  motion  associated  with  the  rigid-body  degree-of-freedom  can  be 
evaluated  as: 


±(9Cg\ 

dt\  89  ) 


e  n-4^  9 


k—0  j—l 


(83) 


For  the  test  article  considered  in  this  manuscript,  Eq.  (13)  yields: 


d_  fdCg\ 
dt  V  99  ) 


=  Ja 


(84) 


where 

Jge  =  (Jg0  +  J'giGi)  (85) 

in  which  the  mass  moment  of  inertia  of  the  encoder  is  also  included. 


C.  Four-Bar  Linkage 


The  Lagrangian  of  the  four-bar  linkage  Ci  is  completely  written  in  terms  of  9  and  9  as  shown  in  Eq.  (44); 
however,  due  to  the  highly  nonlinear  couplings  between  the  angular  quantities,  care  must  be  taken  to  ensure 
that  all  dependencies  are  considered  when  evaluating  the  derivatives  associated  with  Lagrange’s  equation. 
The  contribution  of  the  four-bar  linkage  to  the  equation  of  motion  associated  with  the  rigid  body  degree-of- 
freedom  is  given  by: 


d l  fdCi\  _  dCi^ 
dt  V  d9  )  99 


(86) 


Examination  of  Eqs.  (34)  and  (35),  and  (38)  reveals  that  Si,  S2,  and  C 1  are  functions  of  all  of  the 
angular  position  variables  associated  with  the  four-bar  linkage,  which  can  be  expressed  in  terms  of  the  time 
dependent  motor  angular  position,  i.e.: 


Si  (0(i),a«(0(t)),M0(t)) 

CiWMW)))  (87) 
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For  brevity,  in  the  subsequent  development  the  dependence  of  these  terms  upon  9{t)  is  implied  rather  than 
explicitly  written.  In  order  to  evaluate  Eq.  (86),  it  is  necessary  to  compute  the  time  derivatives  of  these 
functions.  Application  of  the  multi-variable  chain  rule  yields: 

(88) 

(89) 

(90) 

where, 


dSi  _  ~  d6 
dt  1  dt 

dS-2  dO 
dt  2  dt 

dCi  _  ^  dd 
dt  1  dt 


51  = 

52  = 

Ci  = 

However,  from  Eqs.  (34)  and  (35),  it  readily  follows  that: 


r  dSi  8S i  dai  8Si  dcfn  1 

[  80  dai  d0 

'dS2  dS2dat 
89  +  dai  89 
'8Ci  dCidai 

8<t>i  89  \ 
8S2  8<j)i 
8(j>i  89 

89  +  da  88 

Si  = 

s2  = 
Ci  = 


dSi 

~W 

dS2 

de 

dCi 

~d0 


ds, 

£>mSl 


dS2 

dai 

dCi 

dai 


Si 

Si 


dSi 

d(t>i 

dS2 

d<j>i 


(91) 

(92) 

(93) 

(94) 

(95) 

(96) 


where  the  full  expressions  for  the  individual  partial  derivatives  can  be  found  in  the  Appendix.  Hence, 
evaluating  each  of  the  terms  in  Eq.  (86)  yields: 


where 


d 

dt 


dC, 

ae 


Jie0  +  2  <Jie92 


aue- 


Jle  —  (Ju  +  ^2  iS\  +  J3iS2  +  2P1C1S1) 
ale  =  {(PlCl  +  J2iSi)Si  +  J3,S2S2  +  P1S1C1) 


Finally,  Eq.  (86)  becomes: 


d_  fdCA 

dt  ^  88  J 


8  C, 

86 


Jie8  +  &if)82 


(97) 

(98) 

(99) 

(100) 

(101) 


D.  Flexible  Spar 

For  the  Lagrangian  of  the  flexible  spar  in  Eq.  (77),  note  that  all  three  generalized  coordinates,  9,  r/i  and  r/2, 
are  involved.  Hence,  the  flexible  spar  will  contribute  to  three  equations  of  motion.  Evaluation  of  Lagrange’s 
equation  for  the  coordinate  0,  yields: 

dt  ~  ~d¥  =  +  +  Jsed 2  +  crs„d“  +  ksb  (102) 
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where 


Jesg  =  m5|(ci  +  c2r]l  +  c7rfe) 
Jsl  =  mc4S2 
Jsg  =  rhcsS2 

ase  =  mS2S2(c  i  +  c2rij  +  c7r]l) 
kS9  =  2mS29{c2fi\rji  +  c7?)2??  2) 


Evaluation  of  Lagrange’s  equation  for  the  coordinate  rft,  yields: 


where, 


d_  fdcs\ 

dt  V  df]\  ) 


f)C 

— — —  =  jf  9  +  Jg1  i)\  +  as  91  +  ks 

Qrj1  Sr,!  Sm  U  Sr,!  S*71 


Jg  =  mc4S2 

6771 

J!'1  =  mc2 

Ml 

vSrtl  =  rh(c4S2  -  c2S,2?7i) 
=  Eglsc5r 71 


Evaluation  of  Lagrange’s  equation  for  the  coordinate  ry2,  yields: 


d  (dCg 


dt  \di)2 


dC 


37  M  =  JLJ  +  JZV2  +  as92 


dm 


where, 


Js„2  =  ™c8S2 
J!'1  =  mc7 

bV2 

asV2  =  ™(C8^2  -  C7S2T]2) 

—  Esiscior]2 


(103) 

(104) 

(105) 

(106) 
(107) 


(108) 


(109) 

(110) 
(111) 
(112) 


(113) 


(114) 

(115) 

(116) 
(117) 


E.  Complete  Equations  of  Motion 

Finally,  the  complete  equations  of  motion  can  be  written  as: 


Va  +  Jg,+Jla(8)  +  J°.M  -W)  J%(0) 

T 

ah(9)+ase(9,r 71,77:2) 

Jes  (0)  JT  0 

9i 

+ 

<%!  (Mi) 

J9S  (0)  0  JT 

L  *?72  v  J 

jl2_ 

<%2  (M2) 

Ks9(0>Ml.»7l>»?2,J72) 

Qme{ea)  +  Qie(9)  +  Qde{9, 9, 91, 91, 92, m) 

«S„!  (7?l) 

= 

Qsvi  (m)  +  Qdvi  (0, 0, 91,91,92, 92) 

Ks„2M 

Qs^im)  +  Qdg2(9, 9, 91,91,92, 92) 

(118) 


which  can  be  integrated  using  numerical  methods.  The  explicit  dependence  of  the  inertial,  Coriolis,  and 
generalized  force  terms  upon  the  state  variables  is  shown  to  highlight  the  extensive  coupling  present  in  such 
a  system.  The  forcing  term  that  drives  these  equations  is  the  armature  voltage  e0.  Thus,  these  equations 
can  be  used  to  predict  how  the  system  will  behave  when  a  constant  or  time  varying  voltage  is  applied  to  the 
motor  armature.  These  equations  provide  a  powerful  analysis  tool  for  MAV  system  designers.  They  could 
also  be  used  as  a  model  to  support  the  design  of  mechanism  control  laws  that  have  the  potential  to  improve 
the  behavior  of  wing  beat  motion  over  what  could  be  achieved  by  applying  a  constant  voltage  to  the  motor. 

There  are  a  number  of  generalized  forces  that  act  upon  the  system  such  as  a  velocity  dependent  frictional/back- 
EMF  term  associated  with  the  9  coordinate,  which  can  be  written  as  Qie  =  Bg9 ,  the  damping  within  the 
flexible  spar  that  is  proportional  to  the  modal  coordinate  rates  QSvi  =  ^Ci^i9i,  QsV2  =  2£2w2?72,  and  the 
aerodynamic  drag  of  the  spar  that  influences  each  degree  of  freedom  via  Qde,Qdvl ,  and  Qdv2-  Methods  for 
estimating  each  of  these  generalized  forces  will  be  discussed  in  the  next  two  sections. 
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V.  Estimation  of  Generalized  Forces 


The  analytical  method  presented  thus  far  yields  the  form  of  the  governing  dynamical  equations.  There 
are  several  constants  that  are  device  dependent  and  can  only  be  determined  experimentally.  With  respect 
to  the  DC  motor,  the  motor  torque  constant  Kt,  back-EMF  constant  Kb,  and  armature  resistance  Ra  must 
be  determined  experimentally  or  be  provided  by  the  motor  manufacturer.  The  friction  characteristics  of  the 
gear  train  and  linkage  as  well  as  the  modal  damping  of  the  spar  structural  modes  must  also  be  determined 
from  experiments.  A  process  for  the  estimating  each  of  these  terms  is  be  presented  below. 

A.  Motor  Drive  Torque 

The  generalized  force,  Qms  >  is  the  motor  drive  torque,  which  is  given  by: 

Qme  =  (119) 

The  motor  torque  constant  Kt  and  armature  resistance  Ra  must  be  measured,  or  supplied  by  the  motor 
manufacturer.  A  description  of  how  to  estimate  the  motor  torque  constant  from  experiments  is  provided 
below.  The  armature  voltage  ea  forms  the  external  input  into  the  equations  of  motion. 


B.  Generalized  Frictional  Force 

The  overall  generalized  friction  force  Qie  is  a  lumped  parameter  that  includes  the  effects  of  back-EMF  and 
viscous  friction  in  the  gear  train,  linkage  and  the  motor  itself. 


where 


Be 


Qh 


=  BgO 

(120) 

I<TKn  +b) 

Rfri  ) 

(121) 

A  description  of  how  the  terms  in  Eq.  120  were  measured  follows. 


1.  Determination  of  Motor  Characteristics 

The  important  physical  properties  of  most  small  motors  being  used  in  flapping  wing  MAV  prototypes  are 
not  well  characterized  by  their  manufacturer.  This  is  because  many  of  these  motors  are  used  to  power  toys 
or  pager  motors,  and  such  detailed  design  information  is  not  required  for  those  applications.  A  method  for 
extracting  estimates  of  the  required  motor  parameters  is  presented. 

First,  consider  the  problem  of  extracting  of  the  back-EMF  constant  of  the  motor.  By  substituting  Eq.  3 
into  Eq.  5,  one  can  solve  for  Kb'. 

Kb  =  ~l\.ea  —  Raid)  (122) 

0 

Most  of  the  parameters  on  the  right  hand  side  of  the  equation  are  measurable  in  the  laboratory  using  a 
digital  multimeter.  The  measurement  of  9  can  be  performed  using  a  strobe  light,  an  optical  encoder,  or  a 
visual  image  correlation  system.  Alternatively,  Eq.  3  can  be  used  directly  if  an  active  drive  motor  is  coupled 
to  an  optical  encoder  and  to  the  test  motor  via  an  idler  gear.  The  output  voltage  e&  generated  by  the  test 
motor  can  be  measured,  while  the  period  of  the  pulse  train  generated  by  the  optical  encoder  can  be  used 
to  determine  the  motor  speed.  Each  of  the  methods  described  above  have  been  used  by  the  authors  and 
have  been  found  to  produce  nearly  identical  estimates  of  the  back-EMF  constant  on  a  7mm  pager  motor 
manufactured  by  Super  Slicks.  The  motor  used  in  the  test  article  is  known  in  hobbyist  circles  as  a  Red  Super 
Slicks  motor  and  has  an  armature  resistance  of  2.3  Ohms. 

The  estimation  of  the  torque  constant  Kt  followed  the  method  outlined  by  the  Didel  company®.  A  pinion 
was  mounted  on  the  motor  shaft.  An  arm  with  a  splined  hole  near  one  end  was  fitted  over  the  pinion.  A 
point  on  the  opposite  end  of  the  arm  was  oriented  perpendicular  to  the  longitudinal  axis  of  the  arm.  The 
motor  was  mounted  securely  to  a  laboratory  bench  and  the  point  of  the  arm  was  placed  on  a  scale  with  a  0.1 
g  resolution.  Multiplying  the  moment  arm  by  the  force  measurement  yielded  the  stall  torque.  The  torque 

ehttp:/ /www.  didel.com/microkit/moteurs/Motors.html 
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constant  was  then  be  determined  from  Eq.  2.  The  resulting  mean  estimate  for  the  torque  constant  was  6774 
dyne-cm/amp  with  a  standard  deviation  of  619  dyne-cm/amp.  The  high  standard  deviation  was  attributed 
to  the  low  resolution  of  the  scale  used  in  the  test. 

2.  Determination  of  Mechanism  Friction 

An  initial  estimate  of  the  friction  characteristics  of  the  mechanism  was  obtained  by  measuring  the  steady- 
state  speed  of  the  motor,  gear  train,  and  linkage  and  the  current  flowing  through  the  motor.  Note  that  the 
wing  spar  was  not  attached  to  the  mechanism  during  this  test.  A  viscous  friction  model  was  used  to  describe 
the  effect  of  all  sources  of  mechanical  friction  within  the  subsystem,  viz., 

Tf  =  -bO  (123) 

With  knowledge  of  the  motor  torque  constant,  back  EMF  constant,  steady  state  current,  and  steady  state 
speed,  the  viscous  friction  coefficient  was  estimated  from: 

b=-KTia^KrK,  (m) 

@ss  R 

where  9SS  represents  the  steady  state  velocity  of  the  motor.  In  the  work  presented  here,  the  initial  estimate 
of  the  viscous  friction  coefficient  was  obtained  from  experimental  data  from  a  mechanism-only  test.  The 
final  adjustment  of  the  coefficient  was  performed  by  comparing  simulation  time  histories  of  crank  speed  to 
experimental  measurements  of  crank  speed  using  an  optical  encoder.  It  was  found  that  the  viscous  friction 
coefficient  that  caused  the  spar  stroke  periods  of  the  simulation  to  match  those  observed  in  the  experiment 
was  b  =  0.525  dyne-cm/ (rad/sec).  It  should  be  noted  that  this  parameter,  is  very  sensitive  to  the  method 
of  construction  of  the  overall  device.  For  example,  friction  produced  by  the  linkage  rivets  or  pins  is  related 
to  the  skill  of  the  builder  who  must  carefully  deform  the  heads  of  the  rivets  to  ensure  planar  motion  of  the 
linkage  without  producing  excessive  friction. 

C.  Determination  of  Modal  Characteristics  of  the  Spar 

An  estimate  of  the  modal  damping  of  a  spar  or  wing  structure  is  not  normally  available.  A  test  was  conducted 
to  estimate  the  modal  damping  of  the  first  bending  mode  the  spar  used  in  the  experiment.  The  spar  under 
consideration  was  a  carbon  fiber  rod  that  the  manufacturer  markets  as  an  “0.037  inch  diameter”  rod.  The 
rod  was  cut  to  13  cm  in  length.  An  infrared  non-contact  proximity  sensor,  consisting  of  an  LED  and  photo 
transistor,  was  used  measure  the  deflection  of  a  point  on  the  spar,  while  the  root  was  clamped.  The  spar  was 
plucked  and  the  response  was  recorded  at  a  20  KHz  sampling  rate.  The  damping  was  estimated  by  using 
the  well  known  logarithmic  decrement  technique.  The  undamped  natural  frequency  was  also  estimated  from 
the  response.  The  magnitude  of  the  initial  tip  displacement  was  purposely  small,  i.e.  <  1  cm,  in  order 
to  minimize  the  impact  of  aerodynamic  drag  upon  the  damping  estimate.  The  measurement  of  the  spar 
displacement  in  this  test  was  limited  to  sensor  voltage  output,  since  the  measurement  of  frequency  and 
damping  is  not  dependent  upon  the  actual  physical  units  of  displacement,  e.g.  cm. 

Some  important  observations  were  made  while  conducting  the  test.  First,  commercially  available  carbon 
fiber  rod  does  not  have  a  circular  cross  section.  The  cross  section  is  best  described  as  an  oval,  rather  than 
a  “0.037  inch  rod”  with  a  circular  cross  section.  The  cross  section  measurements  varied  by  +/-  .002  inches 
(+/-  .005  cm).  This  means  that  “pluck”  tests  must  be  conducted  with  care,  since  an  initial  deformation 
that  is  not  aligned  with  a  principal  axis  will  produce  transverse  vibrations  along  two  axes,  making  it  difficult 
for  a  non-contacting  proximity  sensor  to  accurately  measure  the  time  history  of  the  displacement,  which  is 
required  for  the  application  of  the  log  decrement  method.  The  best  results  were  obtained  when  the  initial 
deformation  was  aligned  with  a  principal  axis.  The  response,  as  measured  by  the  non-contact  proximity 
sensor,  was  band-pass  filtered  with  break  frequencies  2  times  lower  and  2  times  higher  than  the  theoretical 
estimate  of  the  first  bending  mode  frequency  for  a  spar  with  a  circular  cross  section.  The  power  spectral 
density  of  the  filtered  response  was  used  to  estimate  the  actual  first  bending  mode  frequency. 

Classical  beam  theory  estimates  the  natural  frequency  of  the  first  bending  mode  of  a  cantilever  beam  to 
be  given  by  Eq.  52.  In  order  to  estimate  the  frequency  of  the  first  bending  mode  of  the  spar  from  Eq.  52, 
one  must  use  estimates  of  the  following  spar  properties:  Young’s  modulus,  cross  sectional  area  moment  of 
inertia  ,  mass  per  unit  length,  and  length.  The  mass  of  the  spar  was  determined  to  the  nearest  1  mg  on  a 
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precision  scale  and  the  length  of  the  rod  was  measured  to  the  nearest  .25  mm  using  a  digital  caliper.  The 
mass  per  unit  length  of  the  spar  was  determined  to  be  m  =  0.01085  g/cm.  Using  the  above  information 
and  the  manufacturer’s  specifications  for  Es  =  1.227  x  1012  dynes/cm2,  as  well  as  the  assumption  that 
the  cross  section  of  the  rod  is  circular,  the  idealized  first  bending  mode  frequency  can  be  computed  to  be 
uitheor  —  432.88  rad/sec.  Since  the  product  ESIS  had  a  higher  degree  of  uncertainty  than  the  spar  mass  and 
length,  the  measured  first  bending  mode  frequency  obtained  from  the  peak  of  the  power  spectral  density  was 
used  to  estimate  the  uncertain  product.  The  power  spectral  density  of  the  band-pass  filtered  tip  deflection 
time  history  showed  that  the  natural  frequency  of  the  first  bending  mode  was  located  at  422.79  rad/sec. 
For  a  uniform  circular  cross  section  of  an  0.037  inch  (.09398  cm)  round  rod  with  an  area  moment  of  inertia 
about  the  neutral  axis  of  Is  =  3.83  x  10-6cm4  ,  the  value  of  Young’s  modulus  that  matches  the  observed 
natural  frequency  is  given  by  Es  =  1.1704  x  1012  dyne/cm2,  which  is  slightly  lower  than  the  manufacturer’s 
specification.  The  lower  frequency  is  most  likely  due  to  the  departure  from  the  ideal  circular  cross  section 
of  the  rod.  Indeed,  beam  theory  predicts  variations  of  up  to  ±  5%  in  first  bending  mode  frequency  for  a  rod 
with  an  elliptical  cross  section  where  the  major  and  minor  axes  vary  by  ±  .005  cm  from  a  nominal  .09398 
cm  (.037  inch)  diameter  circular  cross  section.  The  damping  ratio  for  the  first  bending  mode  based  upon  the 
analysis  using  the  log  decrement  method  was  found  be  (q  =  .0025;  the  damping  ratio  for  the  second  bending 
mode  was  assumed  to  be  the  same  as  for  the  first. 

D.  Aerodynamic  Drag  Forces 

The  standard  definition  for  the  generalized  forces  used  in  Lagrange’s  equations  is  given  by:29 

^  =  (125) 

where  Qi  is  the  ith  generalized  force,  Fj  is  the  jth  component  of  the  external  force  applied  to  the  system,  qi 
is  the  itli  generalized  coordinate,  and  is  the  jth  component  of  the  partial  derivative  of  the  position  vector 
to  the  point  of  application  of  the  force  with  respect  to  the  itli  generalized  coordinate.  For  the  derivation  of 
the  generalized  forces  acting  on  the  flexible  spar,  the  differential  contribution  to  the  generalized  forces  is  first 
calculated.  Then  the  differential  contributions  to  the  generalized  forces  are  integrated  to  obtain  generalized 
forces  as  a  function  of  the  angular  position  of  the  motor  and  the  flexible  modal  coordinates.  The  derivation 
that  follows  assumes  that  the  only  external  force  applied  to  the  system  is  due  to  the  aerodynamic  drag 
acting  on  the  flexible  spar  and  that  the  first  two  modal  coordinates  accurately  represent  the  deformation. 
The  flexible  spar  is  modeled  as  a  cylinder  and  a  blade  element  model  is  used  to  determine  the  aerodynamic 
drag.  As  standard  when  using  the  blade  element  method,  the  drag  over  the  cylinder  is  modeled  using  a 
rectangle  with  dimensions  equal  to  the  cylinder’s  length  and  diameter.  Thus,  the  replacement  area  lies  in 
the  Si  —  S3  plane.  Blade  element  theory  states:30 

dFD  =  ^PCdV^  dA  (126) 

where  is  the  magnitude  of  the  velocity  of  a  point  on  the  representative  area  once  it  has  been  resolved  into 
components  that  lie  in  the  plane  of  the  blade  element.  dA  is  the  differential  area  over  which  the  dynamic 
pressure  acts.  Thus  dA  =  dru  dzs.  It  is  convenient  to  describe  the  aerodynamic  force  in  a  frame  that  defines 
the  blade  elements.  The  definition  of  each  blade  element  along  the  spar  should  remain  constant  in  the 
coordinate  frames  in  which  they  are  expressed.  Thus  at  each  point  in  the  representative  area  a  frame  can  be 
defined  that  deforms  along  with  the  spar.  This  deformed  frame  is  related  to  the  frame  associated  with  the 
undeformed  motion  of  the  spar  through  a  single  axis  rotation.  The  deformation  of  the  spar  is  described  in 
Eq.  45.  The  angle  that  defines  the  relationship  between  the  s-frame  and  the  deformed  d- frame  is  described 
by  taking  the  spatial  derivative  of  Eq.  45  with  respect  to  the  ru  variable.  Thus, 

OO 

w'(ru,  t)  =  ^2  &i(ru)Vi(t )  (127) 

i= 1 

As  the  present  analysis  assumes  that  only  two  modes  are  retained, 

w'(ru,t)  =  $i  (ru)m(i)  +  (t)  (128) 
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Rocker  arm  pivot 


Figure  5.  Coordinate  frames  associated  with  undeformed  and  deformed  spar. 
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The  drag  force  acts  in  opposition  to  the  direction  of  the  velocity  defined  in  the  plane  of  the  blade  element. 
Thus,  the  differential  force  vector  is  written  in  the  inertial  frame  as: 


=  ^vldA  (-») 

=  \pCdVt x(-Voo)  dA  (130) 

To  define  begin  with  the  definition  of  the  velocity  of  the  point  written  in  terms  of  the  deformed  frame: 

v$r  =  (-  ($i(r-u)?7i(t)  +  $2 (ru)r)2(t))  fa  cos  w'(ru,t)  +  r^sin  w'(ru,  t)  +  ($i(ru)fn(t)  +  ^2(ru)r]2(t))  sin  w’ (ru,t 

+  ^2[ru)r]2{t))  fa  sin  w'{ru,t)  +  r(j> cos  w' (ru ,  t)  +  ($i(ru)r)i(t)  +  ^2(ru)r]2(t))  cos  w'(ru,  t)^j  d2 

(131) 


As  the  blade  element  lies  in  the  d2  —  d3  plane,  is  defined  by  ignoring  the  d\  component  of  v£r.  Defining 
Vooy  as  the  d2  component  of  and  writing  the  resulting  vector  in  the  inertial  frame  yields: 

Voc  =  (—  cos  fa  sin  w'(ru,t)  —  sin  fa  cos  w'(ru,  t))  V^yhi  + 

(—  sin  fa  sin  w' (ru,  t)  +  cos  fa  cos  w'(ru,t))  Vooyh2  (132) 

The  position  vector  to  any  point  on  the  rectangle  representing  the  spar  can  be  written  as: 

fpr  =  (-r>  cos  fa  -  ru  cos  fa  -  sin  fa$irji  -  sin  fa$2 r]2)  hi 

+  (?>  sin  fa  +  ru  sin  fa  +  cos  fa<&iT)i  +  cos  fa^2r]2)  h2  +  zh3  (133) 


From  Eqs.  130  and  133  it  is  shown  that  both  the  differential  force  and  the  position  vector  from  the  origin 
to  the  point  of  application  of  the  differential  force  are  not  explicit  functions  of  the  generalized  coordinate 
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9 ;  however,  recall  that  (f>i  is  a  function  of  9.  Thus  the  chain  rule  for  partial  derivation  must  be  used  in 
the  determination  of  the  differential  contribution  to  the  generalized  force  associated  with  the  motor  angular 
position.  Thus, 


dQde 


dFpx 


<9rpr.x  d  (pi 
d(j>i  39 


+  dFDy 


3rpry  d^i 

3(j)i  39 


(134) 


dQ^=dFDx^ 

=  dFD,^ 


+  dFDy 

+  dFDy 


dr. 


pry 
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(135) 

(136) 


Thus 
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(138) 


Qdyy  ~ 
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dr. 


dm  1  ““ Dy  dm 

Note  that  these  generalized  forces  are  state  dependent  and  therefore,  these  spatial  integrals  must  be 
evaluated  at  each  simulation  time-step  by  quadrature. 


prv 


(139) 


E.  Summary  of  Measured  Parameters 

Table  1  shows  the  estimates  of  the  parameters  of  associated  with  the  drive  train  elements.  The  inertias  of 
the  gears  were  estimated  by  weighing  each  gear  on  a  precision  scale  to  the  nearest  0.1  mg.  The  gears  were 
constructed  in  a  3D  CAD  program  and  mass  moment  of  inertia  properties  were  estimated  using  the  mass 
measurement  and  the  computed  volume  to  obtain  the  density  of  the  material.  The  inertia  of  the  coreless 
motor  armature  was  estimated  by  dissection  and  measurement.  The  inside  and  outside  diameter  of  the 
copper  windings  were  measured,  clipped  off  and  weighed  on  a  scale  to  the  nearest  0.1  mg.  Similarly  the 
spindle  was  measured  and  weighed  separately.  The  armature  assembly  was  drawn  in  a  3D  CAD  program  and 
the  mass  properties  were  estimated.  It  is  important  to  accurately  characterize  the  armature  inertia  because 
there  is  no  mechanical  advantage  at  the  motor  stage  to  mitigate  its  effect. 


Table  1.  Parameters  for  the  motor,  gears  and  encoders 


Descriptions 

Variables 

Values 

Units 

Stage  1  gear  ratio 

G-\ 

9/48 

- 

Stage  2  gear  ratio 

g2 

12/81 

- 

Inertia  of  stage  1  gear 

JGX 

0.032334 

g-cm2 

Inertia  of  stage  2  gear 

!g2 

0.25515 

g-cm2 

Inertia  of  encoder 

Ie 

0.35329 

g-cm2 

Inertia  of  9-tooth  pinion 

CO 

0.00013719 

g-cm2 

Inertia  of  12-tooth  pinion 

Ip  12 

0.0076397 

g-cm2 

Inertia  of  motor  armature 

la 

0.015845 

g-cm2 

Motor  torque  constant 

I\t 

6774 

(dyne-cm) /A 

Motor  back-EMF  constant 

Kb 

0.000764 

V/ (rad/sec) 

Mechanism  friction  coefficient 

b 

0.525 

dyne-cm/  (rad/sec) 

Motor  resistance 

Ra 

2.3 

n 

Table  2  shows  the  estimates  of  the  parameters  of  associated  with  the  drive  train  elements.  These  estimates 
were  obtained  by  weighing  the  linkage  elements  on  a  scale  to  the  nearest  0.1  mg  and  drawing  the  elements 
in  a  3D  CAD  package  that  provides  mass  property  estimates. 
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Table  2.  Parameters  for  the  four-bar  linkage 


Descriptions 

Variables 

Values 

Units 

Line-of-centers  length 

lo 

-1.936 

cm 

Crank  length  (link  1  formed  by  81-tooth  gear) 

h 

0.461 

cm 

Link  2  pin-to-pin  length 

1-2 

2.060 

cm 

Link  3  pin-to-pin  length 

h 

0.670 

cm 

Distance  to  CG  of  link  2  from  connection  point  at  link  1 

lc2 

\h 

cm 

Distance  to  CG  of  link  3  from  connection  point  at  link  2 

lc  3 

0.753 

cm 

Distance  to  root  of  spar  from  rocker  pivot 

ry 

1.735 

cm 

Mass  of  pins  at  linkage  joints 

rrip 

0.025 

g 

Mass  of  link  2 

m2 

0.1072 

g 

Mass  of  link  3 

m3 

0.1513 

g 

Moment  of  inertia  of  link  2  about  its  CG 

h 

0.06165 

g-cm2 

Moment  of  inertia  of  link  3  about  its  CG 

h 

0.04283 

g-cm2 

Table  3  summarizes  the  parameter  estimates  associated  with  the  wing  spar. 


Table  3.  Parameters  for  the  flexible  wing  spar 


Descriptions 

Variables 

Values 

Units 

Spar  Length 

L 

13.0 

cm 

Young’s  Modulus  of  Spar 

Es 

1.1704  x  1012 

dyn  /  cm2 

Diameter  of  spar 

ds 

0.0940 

cm 

Area  moment  of  inertia  of  spar  about  neutral  axis 

Is 

3.8292  x  10"6 

4 

cm 

Mass  per  unit  length  of  spar 

rh 

0.01085 

g/cm 

Length  of  spar  adjoining  rocker  arm 

lr 

0.556 

cm 

Modal  damping  ratio 

p, 

II 

0.025 

- 

VI.  Results 


A.  Experiment 

A  bench-test  experiment  was  performed  using  the  mechanism  shown  in  Figure  1  for  the  purpose  of  providing 
a  point  of  comparison  for  the  simulation  results.  Four  time-based  measurements  were  acquired  from  the 
experiment,  namely,  the  velocity  of  the  first  stage  gear  via  an  optical  encoder,  armature  voltage,  motor 
current,  and  digital  images  from  a  high  speed  camera.  Motor  voltage,  current,  and  encoder  circuit  waveforms 
were  recorded  by  a  data  acquisition  (DAQ)  system. 

An  optical  encoder  with  16  slots  ran  at  the  same  speed  as  the  first  stage  gear,  which  rotated  at  3/16  of 
the  speed  of  the  motor  and  27/4  times  faster  than  the  speed  of  the  crank  gear.  This  arrangement  yielded  a 
resolution  of  108  pulses  per  crank  revolution  and  3.357  pulses  per  motor  revolution.  The  time  delay  between 
encoder  pulses  was  measured  by  a  20  MHz  clock  within  the  data  acquisition  system.  The  pulse  periods  were 
measured  from  falling  edge  to  falling  edge.  The  detection  of  falling  edge  pulses  also  triggered  A/D  sampling 
of  the  instantaneous  armature  voltage  and  supply  current.  A  2000  frame  per  second  camera  system  captured 
overhead  images  of  the  system  against  the  background  of  a  protractor  plate  that  was  centered  on  the  rocker 
arm  pivot.  Protractor  measurements  for  the  spar  tip  could  be  read  to  the  nearest  1/2°  at  low  tip  velocities; 
however,  the  spar  appeared  blurred  over  arcs  of  up  to  3°  at  high  tip  velocities.  Precise  measurements  of 
deflection  were  therefore  only  taken  at  the  stroke  reversal  points. 

The  experiment  proceeded  by  triggering  the  camera  system,  applying  2.81  V  from  a  regulated  power 
supply  to  the  motor  armature,  and  recording  a  4  second  segment  of  the  system  response.  The  polarity  of 
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the  voltage  applied  to  the  armature  resulted  in  a  clockwise  crank  rotation  when  viewed  by  the  cameras, 
which  is  important  to  note  because  the  system  response  is  direction  dependent.  The  4  second  segment  was 
sufficient  for  the  system  to  achieve  a  long  period  of  steady  state  operationf  .  When  synchronizing  the  camera 
clock  to  the  data  acquisition  system  clock,  the  first  frame  that  motion  was  detected  using  the  camera  system 
was  assumed  to  correspond  to  the  first  frame  acquired  by  the  data  acquisition  system.  The  accuracy  of  the 
time-base  synchronization  is,  therefore,  estimated  to  be  within  ±. 0005s. 

B.  Comparison  of  Simulation  Results  to  Experimental  Measurements 

A  segment  of  the  simulation  and  experimental  results  for  the  system  under  steady-state  operation  is  exam¬ 
ined.  Figure  6  shows  a  comparison  of  simulation  time  histories  and  experimental  measurements  that  have 
been  synchronized  based  upon  extremal  spar  tip  displacements.  The  experimental  measurements  that  are 
presented  here  are  typical  of  results  that  were  observed  over  multiple  experimental  trials.  The  middle  two 
displacement  traces  show  the  simulated  time  history  of  the  angular  displacement  of  the  rocker  arm  (blue) 
and  spar  tip  (green)  as  would  be  measured  by  the  background  protractor  on  the  bench  test  platform.  The 
vertical  red  and  blue  lines  indicate  the  times  at  which  the  extremal  spar  tip  displacements  were  measured 
using  the  high  speed  camera  images.  It  can  be  seen  that  the  period  of  the  experimental  results  matches 
the  period  of  the  simulation  results.  This  match  was  achieved  by  fine-tuning  the  friction  coefficient  b  in  the 
simulation  until  the  crank  rotation  periods  coincided,  leaving  all  other  parameters  fixed  at  their  measured 
or  computed  values.  The  measured  values  of  the  spar  tip  deflection  at  the  extremal  points  were  extracted 
from  the  camera  images  and  appear  as  red  and  blue  asterisks  on  the  vertical  lines  associated  with  stroke 
reversals.  Note  that  the  spar  tip  deflections  associated  with  the  blue  asterisks  are  greater  than  those  as¬ 
sociated  with  the  red  asterisks.  This  is  because  the  linkage  produces  temporally  asymmetric  motion  that 
results  in  a  stronger  reversal  on  one  side  than  the  other.  This  effect  can  be  seen  in  the  displacement  plot 
as  well,  where  the  peaks  associated  with  strong  reversals  are  noticeably  sharper  than  those  associated  with 
the  weak  reversals.  The  simulated  tip  deflections  compare  favorably  with  those  measured  by  the  camera 
at  the  extremal  points.  Finally,  the  experimentally  measured  and  simulated  crank  frequencies  can  be  seen 
in  the  top  two  traces.  It  can  be  seen  that  the  characteristics  of  the  time  histories  compare  favorably.  The 
simulation  slightly  under-predicts  the  peak  experimental  crank  frequency  of  14.85  Hz  by  predicting  a  value 
of  13.98  Hz  for  a  difference  of  4.9%.  The  simulation  also  slightly  under-predicts  the  minimum  experimental 
crank  frequency  of  11.03  Hz  by  predicting  a  value  of  11.66  Hz  for  a  difference  of  5.7%. 

Another  measure  of  the  quality  of  the  simulation  is  to  compare  “multi-exposure”  images  taken  of  the 
spar  as  it  traverses  from  one  reversal  point  to  another.  Figure  7  shows  a  comparison  of  enhanced  camera 
images  and  the  simulation  images  for  the  spar  as  it  traverses  from  the  weak  to  the  strong  reversal  point  for 
clockwise  crank  rotation.  Reflections  and  blurred  images  obscure  the  spar  in  some  segments  of  the  motion, 
thus,  the  images  have  been  enhanced  by  the  use  of  a  sketch  overlay  in  order  to  improve  the  visibility  of  the 
spar.  The  sketched  overlay  was  produced  by  carefully  tracing  over  the  centerline  of  each  spar  image  with  a 
solid  black  digital  “pen”  using  image  processing  software.  The  images  shown  in  Figure  7  correspond  to  the 
portion  of  the  time  histories  between  0.802  sec  and  0.841  sec  shown  in  Figure  6.  The  most  striking  feature  of 
the  images  is  the  “dwell”  that  occurs  near  the  60°  mark.  Another  less  significant  dwell  can  be  seen  between 
95°  and  105°.  This  dwell  effect  is  not  readily  discernable  from  time  history  plots  of  state  variables;  however, 
the  effect  is  a  prominent  feature  of  the  motion  that  can  be  seen  when  viewing  video  images  from  the  high 
speed  camera  system.  The  appearance  of  the  dwell  phenomenon  is  a  result  of  the  combined  motion  of  the 
flexible  structure  and  the  rocker  arm.  In  the  dwell  regions,  the  spar  is  moving  backwards  relative  to  the 
rocker  arm  as  a  result  of  dynamic  structural  deformation.  Clearly,  the  presence  of  such  dwell  regions  have 
a  significant  effect  upon  the  local  dynamic  pressure  experienced  by  the  spar  and  the  resulting  aerodynamic 
loads. 

fDefined  as  a  condition  where  the  period  of  crank  rotation  is  constant  from  one  wing  beat  cycle  to  the  next 
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Figure  6.  Comparison  of  Experimental  and  Simulation  Results 
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(b)  2000  FPS  Camera  (Sketch  overlay) 


Figure  7.  Multi-exposure  comparison  on  downstroke  (weak  to  strong) 
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Figure  8.  Multi-exposure  comparison  on  upstroke  (strong  to  weak) 
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Figure  9.  Naive  simulation,  constant  velocity  motor,  rigid  spar 


The  impact  of  the  modeling  error  associated  with  the  kinematic  simulation  on  aerodynamic  loads  can 
be  assessed  by  comparing  the  behavior  of  two  time  varying  integrals  over  the  course  of  a  wingbeat.  The 
magnitude  of  an  arbitrary  aerodynamic  load,  F,  created  by  a  dynamic  pressure  distribution  over  the  length 
of  the  spar  is  given  by: 

F  =  ^pCpds  J  V2[r u,  t)dru  (140) 

If  over  the  course  of  a  wingbeat,  the  terms  premultiplying  the  velocity  squared  integral  are  taken  to  be 
constant,  then  the  differences  between  the  integrals  for  the  kinematic  simulation  and  the  dynamic  simulation 
are  a  measure  of  the  differences  in  instantaneous  aerodynamic  force.  For  the  full  dynamic  simulation, 

L 

V2(ru,  t)dru  =  c\S\d2  +  c2(S292tjI  +  rfr)  +  2c4S2df]i  +  C7(S%b2r]%  +  rfe)  +  2c8S29i)2  (141) 

whereas  for  the  kinematic  simulation, 

Vk{ru,t)dru  =  CiSf02  (142) 

A  plot  that  compares  the  two  integrals  over  the  course  of  a  wingbeat  is  shown  in  Figure  10.  It  can  be 
seen  that  the  flexible  dynamics  significantly  influence  the  value  of  this  integral  and  therefore  the  magnitude 
of  the  aerodynamic  loads.  It  was  found  that  the  cycle-averaged  value  of  the  integral  from  the  dynamic 
simulation  was  3.1659  x  106cm2/sec2,  while  the  value  corresponding  to  the  kinematic  simulation  results  was 
2.5511  x  106cm2/sec2,  which  indicates  that  the  aerodynamic  loads  produced  by  the  flexible  spar  would  be 
about  24%  higher  than  those  predicted  by  a  simple  kinematic  simulation.  It  can  be  concluded  that  despite 
fluctuations  in  the  motor  speed  and  the  presence  of  dwell  regions,  the  flexible  spar  produces  higher  average 
aerodynamic  loads  than  the  rigid  spar.  The  higher  loads  are  the  result  of  an  increase  in  effective  stroke 
amplitude  and  velocity  near  the  stroke  reversal  points. 
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Figure  10.  Behavior  of  velocity  squared  integral,  related  to  magnitude  of  aerodynamic  forces 


An  additional  observation  was  that  the  inclusion  of  the  dynamics  of  the  second  bending  mode,  in  this 
particular  case,  was  not  warranted.  A  short  study  was  conducted  where  all  terms  in  Eq.  118  that  were 
associated  with  the  second  bending  mode  were  set  to  zero.  The  simulation  results  generated  by  the  model 
that  included  the  rigid  body  coordinate  and  the  first  bending  mode  were  visually  indistinguishable  from 
those  shown  in  Figures  6,  7,  8,  and  10.  Furthermore  the  change  in  the  cycle  averaged  value  of  the  velocity 
squared  integral  occurred  in  the  4th  decimal  place,  i.e.  3.1656  x  106cm2/sec2  for  the  single  mode  versus 
3.1659  x  106cm2/sec2  for  the  second  mode. 

This  work  considers  the  dynamic  coupling  between  ornithopter  subsystems  from  the  motor  to  the  wing 
spar.  In  order  to  analyze  a  system  that  includes  wing  ribs  and  membrane  elements,  at  least  two  additional 
degrees-of-freedom  should  be  considered,  viz.,  wing  torsion,  and  out-of-plane  bending.  An  analysis  that 
includes  these  additional  wing  elements  is  beyond  the  scope  of  the  present  paper  and  will  be  considered  in 
future  work. 
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VII.  Conclusions 


From  the  comparison  of  the  simulation  and  experimental  results,  it  is  concluded  that  the  Lagrangian 
analysis  method  captures  the  major  features  and  trends  observed  in  the  experiment  significantly  better  than 
constant  velocity  motor  and  rigid  spar  kinematic  simulations  used  by  previous  investigators.  Furthermore, 
it  is  concluded  that  dynamic  interactions  between  the  flexible  wing  elements,  drive-train,  and  motor,  have  a 
significant  effect  upon  both  the  instantaneous  and  cycle-averaged  aerodynamic  loads  upon  the  spar.  Thus, 
when  designing  machines  to  optimize  wing  motion  or  attempting  to  match  wing  stroke  patterns  observed 
in  natural  flyers,  one  must  consider  coupling  between  drive  components  and  the  wings  rather  than  simply 
prescribing  wing  kinematics.  The  validation  of  the  equations  of  motion  by  comparing  simulation  results 
to  experimental  measurements,  instills  confidence  in  the  accuracy  of  the  model.  It  was  also  found,  that  in 
the  specific  case  considered  in  the  experiment,  the  second  bending  mode  had  little  effect  upon  the  model 
predictions.  It  is  therefore  expected  that  in  many  cases,  it  may  be  possible  to  ignore  second  bending  mode 
effects  and  still  maintain  an  adequate  level  of  accuracy.  Although  the  model  was  validated  on  a  single  test 
bed,  the  modeling  framework  is  applicable  to  a  large  class  of  MAV  ornithopters.  The  equations  of  motion 
presented  in  this  paper  can  be  used  to  predict  how  DC  motor  powered  flapping  wing  micro  air  vehicle  (MAV) 
subsystem  elements  will  behave  when  a  constant  or  time  varying  voltage  is  applied  to  the  motor  armature. 
These  equations  provide  a  powerful  analysis  tool  for  MAV  system  designers  because  they  can  be  used  as  the 
basis  for  a  systematic  method  of  selecting  motors,  drive-trains,  linkage  kinematics,  and  designing  structural 
and  aerodynamic  properties  of  flapping  wing  MAVs.  They  also  have  the  potential  to  support  mechanism 
control  law  designs  that  can  improve  the  behavior  of  wingbeat  motion  over  what  could  be  achieved  by  simply 
applying  a  constant  voltage  to  the  drive  motor. 
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IX.  Appendix 


A.  Differentiation  of  the  Configuration  Dependent  Terms  in  the  Four-Bar  Linkage 


In  order  to  evaluate  Eq.  (101)  the  following  relationships  are  required: 


dS  i 

dai 


dS2 

d(f>i 


dSi  _  ^2licos{(tn-T0) 
dd  l2  sin(cq  —  fa) 

(143) 

^li  cos(a;  -  4n)  sin(4>i  -  Td) 
l2  sin2  (a;  -  (fn) 

(144) 

dSi  _  ^  h  sin(cq  -  T0) 
d(j)i  l2  sin2(a;  -  <j>i) 

(145) 

dS2  cos(cq  -  T0) 

dd  l3  sin(a;  -  (f>i) 

(146) 

0S2  r  h  sm((f>i  -  T9) 

dai  l3  sin2  (a;  —  (fn) 

(147) 

_  r  k  sin(a;  -  T6»)  cos (cq  -  4>i) 
l3  sin2 (a;  -  (j>i) 

(148) 

-^.  =  -Tsm(T9-ai) 

(149) 
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